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VERIFICATION  AND  VALIDATION  OF  TROPOSPHERIC 
MODEL/DATABASE 


1.  INTRODUCTION 

The  study  of  the  effects  of  the  atmosphere  on  microwave  propagation  necessitates 
a  knowledge  of  height  variation  of  refractive  index  in  tropospheric  and  ionospheric 
regions.  Since  the  magnitude  of  the  refractive  index  is  a  fimction  of  such  parameters  as 
geographic  location,  weather,  time  of  the  day,  and  season  of  the  year,  it  becomes  an 
overwhelming  task  to  analyze  atmospheric  propagation  effects  under  parametric 
conditions.  The  troposphere  research  team  has  developed  an  atmospheric  model,  in  which 
representatives  of  average  climatology(  ECMWF  and  HIRAS  )  data  and  meteorology( 
MRF  and  FNL )  data  are  employed,  for  system  engineers  or  researchers  in  microwave  and 
optical  wave  propagation  field.  Since  the  system  engineer  in  the  propagation  field  needs  a 
few  propagation  parameter  values  for  system  design,  the  researcher  must  develop  the 
propagation  models  either  to  verify  the  model  or  to  generate  a  new  model. 

This  document  summarizes  the  verification  and  validation  techniques  applied  to 
databases  and  models  developed  under  the  Tropospheric  Research  Project  at  the  Naval 
Research  Laboratory.  During  the  course  of  this  project,  several  climatological  data  sets 
are  processed  and  evaluated  by  using  the  ray  trace  algorithm  and  statistical  data  analysis. 
Two  data  formats  are  derived  for  each  data  source.  The  first  data  format  is  the  multilayer 
17-variable  data.  The  second  data  format,  7-variable  surface  data,  is  derived  from  the 
multilayer  data.  Using  the  7-variable  surface  information,  the  Modified  Exponential 
Model^*’^^  builds  a  multilayer  atmosphere.  The  constructed  atmospheric  data  is  then  used 
with  Millman’s  Stratified  Layer  Method^^^  to  calculate  range  error,  time  delay,  and  angle  of 
arrival  error.  In  addition  to  the  surface  climatology  data,  error  parameter  statistics  are 
delivered  under  this  project.  Specifically,  worldwide  time  delay  and  angle  of  arrival  error 
standard  deviations  are  provided.  In  addition,  these  error  standard  deviations  are  sensitive 
to  azimuth  variation. 

This  paper  describes  the  methods,  algorithms  and  safeguards  used  in  verifying  and 
validating  modified  tropospheric  models  with  climatology  databases  such  as  ECMWF  and 
HIRAS.  This  tropospheric  model  code  may  be  easily  interfaced  with  other  users.  There 
are  two  data  files  for  modified  tropospheric  model  such  as  meteorological  and  reference 
height  data.  These  two  data  files  should  generate  more  accurate  time  delay,  range  error, 
and  angle  of  arrival  error  because  these  data  files  consist  of  meteorological  data  and 
reference  coefficient  height  for  each  grid  rather  than  an  average  constant.  Also,  this 
modified  tropospheric  model  can  be  accessible  to  real  time  weather  data. 


Manuscript  approved  November  17,  1997. 
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2.  DATABASE  VERIFICATION 


2.1  Data  Sources 

Climatological  data  are  obtained  from  numerous  sources  including  the  Air  Force, 
the  Navy  and  the  NOAA.  Surface  data  from  the  following  four  databases  is  derived. 

2. 1. 1.  European  Center  for  Medium  Range  Weather  Forecast  (ECMWF) 

ECMWF  data  was  obtained  from  the  US  Air  Force  Environmental  Technical 
Applications  Center  (ETAC)  at  Scott  Air  Force  Base,  Illinois.  The  data  consists  of 
averaged  monthly  over  the  years  1981-1991  with  2.5  x  2.5  degree  grids.  There  are  14 
layers  of  data  by  geopotential  height  and  pressure  levels  from  20  to  1000  millibars 
(mbars).  The  data  elements  include  latitude,  longitude,  temperature,  dew  point,  air 
density,  and  number  of  observations  used  to  obtain  mean  and  standard  deviation  values. 

2. 1. 2  High  Resolution  Analysis  System  (HIRAS) 

HIRAS  data  was  obtained  from  the  ETAC,  Scott  Air  Force  Base,  Illinois.  The 
data  includes  monthly  and  6-hourly  averaged  climatology  data  over  a  nine-year  period 
with  2.5  X  2.5  degree  grids  around  the  world.  There  are  14  layers  data  by  geopotential 
height  and  pre-defined  pressure  levels  from  10  to  1000  mbars.  Other  data  elements 
include  temperature,  relative  humidity,  latitude,  longitude,  date,  and  time.  Associated 
heights  from  the  ECMWF  data  are  used  with  the  HIRAS  data  upon  full  consultations  with 
engineer  and  scientists  who  generated  the  data. 

2. 1. 3  Medium  Rcmge  Forecast  (MRF) 

MRF  data  from  the  NCDC,  NOAA,  includes  6  hourly  diurnal  meteorological  data 
over  the  period  of  1st  January  1991  to  31st  December  1996,  with  2.5  x  2.5  degree  grids 
around  the  world.  There  are  12-layers  from  the  surface  to  the  height  of  50  mbars.  The 
data  elements  include  geopotential  height,  air  temperature,  relative  humidity,  wind  vector 
and  speed,  latitude,  longitude,  date,  and  time. 

2. 7. 4  Final  Analysis  (FNL) 

FNL  Data  are  the  version  of  the  improved  MRF  Data  consisting  of  6  hourly 
diurnal  meteorological  data  over  the  period  of  1st  January  1997  to  present  with  1.0  x  1.0 
degree  resolution  around  the  world.  There  are  13  layers  from  the  Earth’s  surface  to  the 
height  of  20  mbars.  Data  elements  include  geopotential  height,  temperature,  relative 
humidity,  total  cloud  cover,  wind  vector  and  speed,  geopotential  height  and  pressure 
vertical  velocity. 
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2.2  Processing  Multilayer  Data 

ECMWF  data  are  delivered  in  unformatted  ASCII  text.  HIRAS  data  are  delivered 
in  formatted  ASCII  text.  MRF  and  FNL  are  converted  into  ASCII  data  from  binary.  All 
four  data  sources  are  converted  into  the  troposphere  global  17-variable  multilayer  format. 
See  appendix  A  for  a  description  of  the  troposphere  database  format.  ECMWF,  HIRAS, 
and  MRF  data  sets  were  defined  for  each  of  2.5°  x  2.5°  grids;  10,368  total.  FNL  dataset 
is  the  improved  MRF  data.  The  resolution  of  FNL  data  are  1°  x  1°  resolution  as  opposed 
to  the  2.5°  X  2.5°  resolution  of  the  other  data  sources.  These  data  are  checked  to  ensure 
valid  ranges.  Minimum  and  maximum  field  values  are  calculated  and  checked  for 
reasonableness. 

Occasionally,  a  missing  value  is  appeared  in  the  raw  data  for  a  parameter  such  as 
relative  humidity.  Missing  data  are  indicated  consistently  by  a  zero  value  where  a  zero 
value  was  unreasonable.  For  instance,  the  Table  1  presents  a  sample  dataset  with  a  zero 
relative  humidity  value  at  850  mbars  pressure  level.  It  is  expected  that  relative  humidity 
would  decrease  with  increasing  height  and  with  decreasing  pressure.  A  linear  fit  using  the 
two  nearest  temperatures  and  pressures  would  be  used  to  interpolate  the  relative  humidity 
at  the  850  mbars  pressure  level  and  the  surface  pressure  level. 

Table  1 .  Sample  Data  Anomaly  for  ECMWF,  HIRAS  and  MRF 


Pressure 

(mbars) 

Temperature 

(°k) 

Relative 

Humidity 

(%) 

Geopotential 

Height 

(m) 

20 

221 

0 

26354 

30 

216 

0 

23734 

50 

207 

0 

20603 

70 

199 

0 

18624 

100 

198 

0 

16565 

150 

206 

0 

14181 

200 

218 

0 

12387 

250 

230 

0 

10918 

300 

239 

37.6 

9666 

400 

255 

42.8 

7577 

500 

267 

47.3 

5870 

700 

282 

52.6 

3160 

850 

293 

0 

1520 

1000 

301 

63.4 

104 

Although  pressure  layer  data  are  not  expected  to  be  strictly  increasing  or  decreasing,  all 
10,368  grids  data  are  read  and  anomalous  records  written  to  a  file  for  review  and 
assessment.  Abnormal  data  in  temperature  and  height  are  considered  to  be  strictly 
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increasing/decreasing.  Anomalous  relative  humidities  are  out  of  range  and  negative  height 
or  “zero”  data  are  not  reasonable  for  each  level.  Each  of  the  anomalous  records  is 
reviewed  and  corrected  through  interpolation  process.  If  the  value  was  determined  to  be 
missing,  a  simple  linear  interpolation  corrects  the  missing  value.  MRF  and  FNL  data  sets 
are  in  polar  coordinates.  In  order  to  assign  our  coordinates,  grid  system  is  used  to  unpack 
and  convert  the  binary  data.  Additional  code  has  been  developed  to  use  linear 
interpolation  to  fill  the  missing  grids  that  are  resulted  from  polar-to-cartesian  coordinate 
conversion. 

2.3  Validation  and  Verification  of  Multilayer  Data 

The  17  variable  multilayer  data  are  contain  several  occurrences  of  negative  height 
at  1000  mbars  pressure  level.  This  is  due  to  the  fact  that  surface  pressure  is  less  than  1000 
mbars  at  certain  points.  In  order  to  use  the  data  for  raytracing,  pressure  is  interpolated  for 
zero  height  where  a  negative  height  occurred.  In  addition,  data  are  checked  for  zero 
temperature  and  zero  relative  humidity  for  700  mbars,  800  mbars,  and  1000  mbars  and 
interpolated  these  values  whenever  necessary.  The  program  then  calculates  refi'activity 
and  dry/wet  heights  for  the  new  readings.  As  a  final  effort,  the  global  database  is  checked 
to  ensure  that  14  multilayers  occurred  for  each  reported  grid  and  no  negative  heights,  zero 
temperatures  or  zero  relative  humidities  occurred  above  400  mbars.  Worldwide  contour 
maps  are  plotted  for  pressure,  reffactivity,  temperature  and  relative  humidity  for  different 
data  sources  and  visually  inspected  for  patterns  as  well  as  reasonableness.  Figure  la  and 
lb  show  worldwide  contour  maps  for  ECMWF  and  MRF  data  source  of  refi'activity. 

2.4.  Surface  and  Reference  Coefficient  Data 

The  processed  17  variable  multilayer  data  need  lots  of  memory  for  storage. 
ECMWF  data  need  about  6  Mbytes  for  imcompressed  data.  HIRAS  data  need  4  times 
more  than  ECMWF  data  about  24  Mbytes  for  uncompressed  data.  MRF  data  need  about 
730  Mbytes  for  one  year  uncompressed  data.  FNL  data  need  about  1.2  Gbytes  for  one 
month  uncompressed  data.  It  is  inappropriate  to  use  the  processed  17  variable  multilayer 
data.  Therefore  the  seven  variable  surface  data  is  generated  from  the  17  variable 
multilayer  data.  Surface  data  is  defined  as  data  where  height  equals  zero.  Surface  data  is 
provided.  The  data  is  a  result  of  setting  height  equal  to  zero  and  interpolating  the 
remaining  elements.  The  interpolation  process  uses  a  second  order  polynomial  fit  which 
seems  to  produce  the  closest  approximation  to  the  multilayer  data.  The  surface  database 
is  organized  in  a  series  of  2.5®  x  2.5®  grids.  The  grids  are  numbered  beginning  with  grid 
number  1  at  0  degree  latitude  and  0  degree  longitude(  Greenwich,  England)  and 
subsequent  grids  follow  the  equator  in  an  easterly  direction  fi'om  0  degree  longitude  to 
357.5  degree  longitude.  The  second  row  begins  at  2.5  degree  latitude  and  0  degree 
longitude  and  encircles  the  earth  in  an  easterly  direction.  This  continues  to  87.5  degree 
latitude  and  357.5  degree  longitude  for  northern  hemisphere.  There  is  no  data  for  the  ring 
at  exactly  90  degree  north.  The  next  series  of  grids  begin  at  -2.5  degree  latitude  and  0 
degree  longitude  and  terminate  at  -90.0  degree  latitude  and  357.5  longitude  for  the 
southern  hemisphere.  Figure  2  shows  the  database  organization  of  each  data  (  ECMWF, 
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Figure  1b-  MRF  Multilayer  Refractivity(N)  at  lOOOmb  Contour  Map  for  August  9  at  1800 


325  365 


HIRAS,  and  MRF).  For  example,  the  grid  number  1  contains  areas  that  are  covered  in  the 
area  of  greater  than  or  equal  to  0  degree  latitude  and  less  than  2.5  degree  latitude  and  are 
covered  in  the  area  greater  than  or  equal  to  0  degree  longitude  and  less  than  2.5  degree 
longitude.  The  grid  number  2  continues  areas  that  are  covered  in  the  area  of  greater  than 
or  equal  to  0  degree  latitude  and  less  than  2.5  degree  latitude  and  all  longitudes  that  are 
greater  than  or  equal  to  2.5  degree  longitude  and  less  than  5.0  degree  longitude.  This 
logic  is  followed  for  the  remaining  grids.  Table  2  contains  a  sample  of  surface  data. 
Pressure  (Press),  reffactivity  (N),  temperature  (T),  and  relative  humidity  (RH)  are 
interpolated  for  each  day/time  period  for  each  of  the  four  data  sources.  Figure  3  through 
5  show  surface  reffactivity  contour  maps  for  ECMWF,  HIRAS,  and  MRF  dataset. 


Table  2.  Surface  Data  format  for  each  data  file  (ECMWF,  HIRAS,  MRF,  and  FNL) 


GRID  NUMBER 

1  through  10,368  Note:  latitude  90  deg  is  not  report. 

LATITUDE  (  degree) 

Latitude  in  positive  degrees  for  North,  negative  degrees 
for  South 

LONGITUDE  (degree) 

Longitude  in  degrees  East,  0  to  360  degree 

SURFACE 

PRESSURE  Imbars) 

Pressure  in  millibars 

Mean  total  Refractivity  in  N  units  ( Ndry  +  Nevrt) 

SURFACE 

TEMPERATURE  (K) 

Mean  Temperature  in  Kelvins 

SURFACE 

RELATIVE  HUMIDITY  (%) 

Mean  percentage  relative  humidity 

After  the  surface  data  is  used  in  the  Modified  Exponential  Model  to  build  a  local 
atmosphere,  Millman’s  Stratified  Layer  method  of  raytracing  uses  the  constructed  data  to 
produce  the  error  parameters  of  interest  such  as  range  error,  time  delay  and  angle  of 
arrival  error. 

The  reference  height  coefficient  data  consists  of  hourly/daily/monthly/  height 
coefficients,  undulation  value,  time  ,range,  and  angle  of  arrival  error  standard  deviation  of 
the  tropospheric  propagation  model  for  each  grid.  These  coefficients  are  used  to  build  the 
tropospheric  reffactivity  to  about  20  mbars  pressure  level  ffom  the  surface  data.  The 
undulation  of  geoid  is  a  distance  between  the  geoid  and  the  mathematical  reference 
ellipsoid  as  measured  along  the  ellipsoidal  normal.  This  distance  is  considered  positive 
outside,  or  negative  inside,  the  reference  ellipsoid.  Table  3  shows  the  reference  coefficient 
to  each  month. 
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Figure  3-  ECMWF 
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Figure  4-  MRF  Surface  Refractivity(N)  Contour  Map  for  August  15  at  0600 
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Figure  5-  HIRAS  Surface  Refractivity(N)  Contour  Map  for  August  at  0600 


Table  3.  The  format  of  reference  coefficient  datable 


#  of  Grid 

Month 

Undulation 

Standard 

Deviation 

Time/Ranae  Error 

Standard 
Deviation 
Anale  Error 

1 

99.99  Km 

±  20.0  m 

3.03 

0.23 

10,368 

89.00  km 

±  13.6  m 

4.32 

0.06 

3.  MODEL  VERIFICATION 
3.1  Verification  Strategy 

Millman’s  Stratified  Layer  method  of  raytracing  divides  the  troposphere  into 
discrete  layers  of  reffactivity  ni.  The  assigned  value  may  come  from  actual  meteorological 
measurements  (our  multilayer  data)  or  from  modeled  atmospheres.  Three  models  are 
delivered  under  this  project.  The  Modified  Exponential,  Hopfield^'’’^^  and  ESS  Models^®’^ 
are  compared  with  computed  error  parameters  and  actual  measured  multilayer  data. 
Error  parameters  of  interest  are  time  delay,  range  error  and  angle  of  arrival  error.  Several 
areas-of-interest  (AOIs)  were  selected  to  examine  the  performance  results  of  each  model. 
Approximately  130  AOIs  were  selected  for  unique  geographic  and  climatic  conditions. 

3.2.  Modified  Exponential  Model 

Visual  inspection  of  the  multilayer  data  indicates  that  refractivity  is  a  decreasing 
exponential  function  of  height.  In  fact  any  local  refractivity  N(h)  can  be  expressed  as 

N(h)=Nsexp(-h/H) 

where  H  is  a  scale  (or  reference)  height  appropriate  to  the  value  of  N  at  zero  height.  Ns. 
Considering  a  scale  height  here,  it  is  simply  the  height  at  which  the  value  of  N  (h)  is  1/e  of 
the  surface  refractivity  Ns  from  the  above  equation ,  at  which  the  height  h  is  equal  to  the 
scale  height  H.  The  wet  refractivity,  Nw,  is  below  1.0  N-unit  in  comparison  with  the  dry 
refractivity  Na  with  100  N-unit  in  the  neighborhood  of  the  selected  value  of  a  reference 
height,  H.  The  ratio  of  dry  to  wet  refractivity  is  approximately  100  at  the  reference  height 
H  where  the  height  from  the  surface  is  1/3  of  the  total  tropospheric  region.  The  refractive 
effect  of  bending  and  time  delay  beyond  this  layer  (reference  height)  will  be  very  limited 
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since  temperature  and  humidity  does  not  change  drastically  to  affect  refractive  bending  in 
those  extended  areas  such  as  tropopause,  stratosphere,  stratopause,  free  space  and 
ionosphere  above  1  GHz.  This  coincides  with  the  fact  that  the  most  bending  and  refractive 
phenomena  occur  within  this  region  (beneath  the  reference  height)  from  the  surface  of  the 
Earth.  This  implies  that  the  tropospheric  effects  on  the  ray  bending  can  be  approximated 
with  the  reference  height  without  loss  of  any  significant  physical  principle. 

In  an  inversion  layer  temperature  increases  with  altitude,  and  such  a  layer  is  highly 
stable.  All  vertical  motions  are  strongly  inhibited  in  an  inversion  layer,  and  pollution 
emitted  below  the  layer  tends  to  be  confined  below  it.  Also,  if  a  source  of  water  vapor 
exists  below  an  inversion  layer,  it  tends  to  be  confined  below  the  layer,  with  the  result  that 
large  decreases  in  index  of  refraction  may  be  encountered  in  the  upward  passage  through 
an  inversion  layer.  Thus,  the  occurrence  of  inversion  layers  has  an  important  effect  on 
earth-space  communication  paths  in  the  low  elevation  angle.  The  decrease  or  change  of 
the  water  vapor  pressure(  e  )  with  height  is  generally  variable  but  may  be  approximately 
exponential.  Note  also  that  the  delay  caused  by  water  vapor  is  considerably  smaller  than 
that  for  dry  air  above  3  km  to  5  km  from  the  surface,  but  total  water  vapor  content  along  a 
path  is  variable  and  not  predictable  with  high  accuracy  from  the  surface  water  vapor 
pressure  or  density.  Therefore,  water  vapor  is  responsible  for  a  larger  error  or  uncertainty 
in  the  range  calculation  than  in  dry  air  at  lower  atmosphere.  This  kind  of  exponential 
model  is  widely  applicable  and  is  reliable  when  any  reliable  climatological  or 
meteorological  data  on  actual  refractivity  profiles  are  applied.  In  this  report  a  global  2.5^ 
X  2.5®  grid  accuracy  is  used.  This  model  provides  the  accuracy  of  less  than  1%  of  root- 
mean  square  error  (rms)  from  the  climatology  or  meteorology  data  in  comparison  with  the 
accuracy  of  20  %  to  30  %  of  rms  errors  for  the  Hopfield  model  and  other  models. 
Therefore,  this  model  approach  has  been  chosen  in  this  study  as  the  most  reliable  and 
accurate  model  in  comparison  with  other  models  for  various  different  conditions. 

3.3.  Hopfield  Model 

The  actual  atmosphere  is  not  isothermal,  nor  is  its  composition  an  ideal  gas,  and 
the  state  of  water  vapor  in  the  atmosphere  is  quite  different  from  the  state  of  an  ideal  gas. 
So  it  is  more  practical  to  express  the  refractivity  in  the  quadratic  term  induced  by  the  dry 
air  or  the  water  vapor  separately  as  follows: 

N(h)  =  Nd  +  Nw 

where  Nd  is  the  dry  refractivity  and  Nw  the  wet  refractivity  represented  by 

Nd  =  kd(hod  -h)" 

Nw  =  kw(how-h)" 


with  kd  =  l/[hod  -  hs]"*,  kw  =  l/[how  -  hs]'*,  hod  the  dry  height  of  the  order  of  40  km,  and 
how  the  wet  height  of  the  order  of  12  km.  The  hs  is  the  surface  height  and  h  is  the  height 
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of  each  layer.  It  is  noted  that  if  the  reff  activity  as  a  fianction  of  height  is  represented  by  an 
exponential,  it  is  not  integrable  in  closed  form,  where  if  it  has  the  form  as  in  dry  and  wet 
reffactivity  equations  ,  it  is  integrable.  The  representation  of  the  dry  and  wet  terms  of 
refractivity  by  quadratic  equations  as  dry  and  wet  reffactivity  equations  give  a  good 
agreement  with  range  error  and  Doppler  data  above  6®  elevation  angles.  Note  here  that  if 
monthly  or  weekly  averages  of  the  reffactivity  are  used  and  the  accuracy  is  not  so 
important,  this  model  is  useful  and  practical. 


3.4.  ESS  Model 


The  ESS  Model  is  a  second  order  polynomial  equation  with  varying  zenith  angle 
(za)  from  ground  station  to  satellite.  Negative  angles  are  converted  to  0  degree  and 
below  the  horizon  measurements  (  >  90  degrees)  are  converted  to  90  degrees.  The 
equation  for  a  range  delay  is: 


D 

RANGEDEL  =  - 

A  +  B*cos(za)  +  C  *  cos^  ( za  ) 

where  D  is  equal  to  0.3048006. 


A  =  0.003589 

B  =  0.087605 

C  =  0.19696793 

za  <5° 

A  =  0.002129 

B  =  0.12158 

o 

II 

U 

5‘’<za<30' 

A  =  0.030000 

B  =  0.080000 

o 

II 

CJ 

za<  30° 

4.  MODEL  VALIDATION 
4.1.  Selected  Area-of-Interest 

Raytrace  error  parameters  are  dependent  on  the  time  and  the  place.  The  130  area- 
of  -interest  are  selected  as  the  representative  geographic  and  climatic  global  conditions  for 
the  completion  of  model  effectiveness.  Table  4  contains  a  list  of  these  areas  of  interest 
with  climatic  and  topographic  descriptions.  Figures  6  and  7  show  time  delay  and  angle 
error  calculated  ^vith  the  Modified  Exponential,  Hopfield  and  ESS  models  using  surface 
data  plotted  against  the  actual  ECMWF  data  on  August  and  February  data  at  Tehran,  Iran. 
The  results  show  that  the  Modified  Exponential  Model  is  closest  to  the  measured  data, 
and  has  been  consistent  regardless  of  area-of-interest,  month,  or  data  source. 
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Table  4  .  Selected  Areas-of-Interest 


Climate 

1 

LU 

1 

TEMPERATE 

TEMPERATE 

1 

LU 

CL 

BOREAL 

SUBTROPICAL 

TEMPERATE 

TROPICAL 

I 

LU 

CL 

TROPICAL 

1 

LU 

CL 

TROPICAL 

BOREAL 

TROPICAL 

TROPICAL 

SUBTROPICAL 

i 

CL 

SUBTROPICAL 

TROPICAL 

TROPICAL  1 

TEMPERATE 

TEMPERATE 

POLAR 

BOREAL 

TEMPERATE 

1 

LU 

CL 

TEMPERATE 

Continent 

N.  America 

N.  America 

N.  America 

N.  America 

WATER/Pac 

N.  America 

Europe 

WATER/Med 

WATER/Med 

o 

LU 

1 

WATER/Pac 

N.  America 

N.  America 

N.  America 

S.  America 

Africa 

Africa 

Australia 

Asia 

Asia 

Asia 

<D 

CL 

O 

1- 

D 

LU 

Asia 

N.  America 

N.  America 

N.  America 

N.  America 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

HI  LON 

o 

d 

o 

CO 

260.0 

o 

d 

o 

CO 

o 

d 

CO 

CM 

190.0 

285.0 

CD 

LO 

o 

d 

o 

d 

90.0 

155.0 

155.0 

o 

d 

CO 

290.0 

325.0 

35.0 

O 

d 

I 

I 

I 

LO 

cvi 

90.0 

o 

d 

00 

195.0 

292.5 

292.5 

o 

d 

CD 

CM 

LO  LON 

o 

d 

CO 

CM 

225.0 

285.0 

245.0 

165.0 

275.0 

345.0 

30.0 

345.0 

60.0 

a 

d 

135.0 

o 

d 

CO 

CM 

260.0 

285.0 

75.0 

105.0 

85.0 

CD 

d 

CO 

o 

d 

CD 

CD 

d 

287.5 

O 

d 

O) 

CM 

287.5 

Hi  LAT 

o 

d 

LO 

55.0 

50.0 

45.0 

60.0 

32.5 

o 

d 

45.0 

50.0 

5.0 

50.0 

22.5 

CD 

d 

CO 

22.5 

o 

d 

CM 

1 

o 

d 

CO 

O 

d 

1 

20.0 

o 

d 

CM 

47.5 

60.0 

o 

d 

00 

o 

d 

CD 

47.5 

42.5 

42.5 

LOLAT 

LO 

cvi 

CM 

22.5 

o 

d 

'M' 

o 

d 

CO 

45.0 

o 

CO 

35.0 

1 

27.5 

1 

22.5 

5.0 

47.5 

I 

1 

1 

1 

37.5 

o 

d 

o 

d 

CD 

45.0 

42.5 

o 

d 

CD 

d 

Description 

Eastern  US 

Western  US 

Northeast  US 

Midwest  US 

Bering  Sea 

Southeast  US  Region  3 

Europe 

Persian  Gulf 

Mediterranean 

Mid-Indian  Ocean 

East  China  Sea  &  Sea  of  Japan 

Northwest  Pacific 

Canada  Belt 

Central  America 

Amazon  Forest 

South  Africa 

Sahara  Desert 

Australia  Continent 

Southeast  Asia  Region  1 

Southeast  Asia  Region  2 

Himalayas 

Eurasia  Belt 

Siberia 

New  Alaska 

State  of  Maine 

Boston,  Massachusetts 

New  York,  New  York 

AOI 

EUS 

WUS 

NEUS 

MWUS 

AK 

SEUS3 

BJR 

e 

MED 

MIO 

Hi 

NWP 

CAN 

CAM 

AMFOR 

SAF 

SAH 

snv 

SEAS1 

CM 

LU 

CO 

GOBI 

EURAS 

SIB 

NAK 

MAINE 

BOSTN 

OAN 

B 

CM 

CO 

B 

LO 

CO 

B 

CO 

G) 

o 

- 

CM 

CO 

CO 

CD 

N- 

CD 

o> 

20 

CM 

1  22 

23 

I  24 

25 

26 

27 

15 


28  OCNCY  Ocean  City,  Maryland _ 37.5  40.0  285.0  287,5 _ N.  America  TEMPERATE 

29  VABCH  Virginia  Beach,  Virginia _ 35.0  37.5  282.5  285.0 _ N.  America  SUBTROPICAL 

30  MYRBC  Myrtle  Beach,  South  Carolina _  32.5  35,0  280,0  282.5 _ N.  America  SUBTROPICAL 

3 1  JAX  Jacksonville,  Florida _ 30.0  32.5  277.5  280.0 _ N.  America  SUBTROPICAL 

32  MIA  Miami.  Florida _ 25.0  27.5  277.5  280.0  IN.  America  |  SUBTROPICAL 
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6 1  HAAM  Himalayas _ 32.5  42.5  245.0  257.5  >3000  Asia _ TEMPERATE 

62  NDSD  North  &  South  Dakota _ 42.5  50.0  255.0  265.0 _ N.  America  TEMPERATE 

63  WCUS  West  Cost  US:  Washington  &  Oregon  42.5  50,0  235.0  240.0 _ N.  America  TEMPERATE 

64  DC _ Washington.  D.C. _ 35.0  40.0  280.0  285,0 _ N.  America  TEMPERATE 

65  NWAFR  West  Africa:  Northwest  Coast  2.5  17.5  340.0  355.0  Africa  TROPICAL 
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Elevation  Angle  Tehran,  Iran  (August,  ECMWF  Data) 


Elevation  Angle  (deg) 


4.2.  Error  Parameter  Statistics 


Particular  interest  to  this  project  is  the  derivation  of  time  delay,  range  error,  and 
angle  of  arrival  error  for  positive  elevation  angle  and  negative  elevation  angle  from  the 
horizontal  direction.  Also  the  standard  deviations  of  each  error  parameters  are  derived. 
The  calculation  method  of  standard  deviation  is  discussed  below.  In  addition,  azimuth 
sensitivity  is  analyzed. 

4. 2. 1.  Time  Delay,  Range  Error,  and  Angle  of  Arrival  Error  for  Negative  Elevation 

In  a  satellite  communication  sytem,  the  transmitter  or  receiver  antenna  of  the 
ground  stataton  transmits  or  receives  signal  above  90  degree  form  zenith  angle  or  negative 
angle  from  horizental  direction.  The  modified  exponential  tropospheric  model  generates 
time  delay,  range  error,  and  angle  of  arrival  error  for  any  negative  elevation  angle  between 
a  ground  station  and  satellite  station.  Figure  8a,  8b,  and  8c  show  the  ECMWF  worldwide 
contour  maps  of  time  delay,  range  error,  and  angle  of  arrival  error  for  90.5  degree  zenith 
angle  as  an  example. 

4.2.2,  Standard  Deviations 

Percent  standard  deviations  are  calculated  for  time  delay  and  zenith  launch  angle 
for  each  of  the  10,368  grids.  For  each  grid  TDu,  ray  traces  were  run  for  the  target  grid 
and  surrounding  grids  with  the  grid  crossing  feature  turned  off.  Table  5  shows  grid  time 
delays  required  to  calculate  standard  deviation  for  TDl. 

Tables.  Sample  Grid  for  Standard  Deviation  Calculation 


TD2 

TD3 

TD4 

TD5 

TDl 

TD6 

TD7 

TD8 

TD9 

First,  the  time  delay  is  determined  for  each  of  the  surrounding  grids.  Then  the  average 
and  standard  deviation  is  calculated  for  TDli; 


TDavg  =  ITDi  /  m 

m 

TDvai  =  Z  (TDi  -  TDavg)^  /  m  -  1 


Figure  8b  -EC 
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99.0 


Figure  8c  -ECMWF  Angle  of  Arrival  Error  Contour  Map  for  March  ,  Truezenang  =  90.5  deg. 
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TDstdev  ~  V  TDvar 


where  m  is  9  (plus  or  minus  one  surrounding  grid)  for  grids  between  60°  North  and  60° 
South  in  latitude.  For  latitudes  between  60°  and  70°  North/South,  plus  or  minus  two  grids 
from  the  center  grid  covers  25  grids.  For  latitudes  between  70°  and  80°  North/South, 
plus  or  minus  three  grids  from  the  center  of  grid  covers  49  grids.  For  latitudes  between 
80°  and  90°  North/South,  plus  or  minus  four  grids  form  the  center  of  grid  covers  8 1  grids. 
Next,  percent  standard  deviation  is  taken  by  dividing  the  standard  deviation  by  the  original 
time  delay  and  multiplying  by  100. 

%  StDev  =  (TDstdev  /  TDl )  *  100 

These  calculations  are  performed  for  all  grids  for  12  months  using  ECMWF  data  and 
HIRAS  data.  Appendix  D  show  contour  plots  of  percent  standard  deviation  for  time 
delay  using  ECMWF  data. 

4.2.3.  Azimuth  Sensitivity 

The  time  delay  is  calculated  for  every  other  grid  of  the  10,368  grids.  Time  delay 
for  the  case  of  azimuth  independent  is  compared  with  that  of  azimuth  dependent.  The 
range  between  minimum  and  maximum  time  delays  for  the  azimuth  dependent  case  is 
about  45  %  of  those  in  the  azimuth  independent  case.  Contour  maps  for  azimuth 
sensitivity  referencing  data  and  fflRAS  data  are  shown  in  appendix  E  respectively. 

5.  VALroATION  OF  THE  MODIFIED  EXPONENTIAL  MODEL  CODE 

The  modified  exponential  model  code  is  designed  as  a  subroutine  package  which 
can  be  incorporated  into  an  application.  All  code  is  written  in  ANSI-standard  FORTRAN 
77  with  no  machine-dependent  extensions  except  VAX  system  (the  data  files  in  VAX 
system  need  a  special  format  for  direct  access  data  ).  This  program  is  divided  into  two 
major  subroutine.  One  is  I/O  subroutine(getwxfil)  and  the  other  is  the  main  calculation 
routine(raytrace).  These  two  main  subroutines  are  further  divided  into  two  subroutines 
in  order  to  calculate  time  delay,  range  error,  and  angle  of  arrival  error. 


The  tropospheric  model  code  would  invoke  the  user  information  and  user  request 
database.  Then,  it  will  proceed  to  call  main  calculation  routine  with  user  information. 
This  model  program  will  generate  error  or  status  messages  to  validate  each  subroutine  and 
function.  If  an  error  or  a  status  flag  is  returned  from  a  subroutine  or  function,  the 
troposphere  program  will  terminate  or  return  an  error  message  indicating  the  subroutine 
where  it  was  generated. 


25 


5.1.  Main  Program 

The  main  program  initializes  and  defines  variables  needed  in  I/O  subroutine  and 
raytrace  calculation  subroutine.  I/O  subroutine  reads  two  data  files.  One  contains 
weather  data  containing  grid  number,  latitude,  longitude,  surface  reffactivity,  pressure, 
temperature,  and  undulation  of  geoid  parameter.  The  other  contains  reference  height 
coefficient  for  each  2.5°  x  2.5°  degree  grid.  The  user  can  choose  one  option  whether  grid 
cross  or  without  grid  cross  features.  If  the  grid  cross  feature  contains  new  grid 
information  (  surface  refractivity,  pressure,  and  temperature)  of  a  ray  path  crossing,  new 
grid  information  rather  than  previous  grid  information  will  be  used  to  calculate  the 
refractivity  with  its  respective  height.  If  a  ray  path  crosses  more  than  two  grids,  raytrace 
subroutine  will  repeat  above  process  to  calculate  time  delay,  range  error,  and  angle  of 
arrival  error.  Raytrace  subroutine  will  use  the  arguments  passed  from  the  main  routine  to 
generate  time  delay,  range  error,  and  angle  of  arrival  error.  The  mathematical  equations 
used  for  this  routine  are  described  in  section  5.3.  The  user  input  is  as  follows: 

(Q)  Enter  the  grid  cross  feature  (  0  :  without  grid  cross,  1  ;  with  grid  cross) 

(O)  Enter  Source  option  ( 1, 2,  3,  or  4  ) 

1.  ECM  database 

2.  MRF  database 

3.  HIRAS  database 

4.  User  database 

(Q)  Enter  the  location  of  station  ( latitude,  longitude  ) 

(Q)  Enter  month,  day  or  time 

(Q)  Enter  the  station  height 

(Q)  Enter  the  Elevation  angle  (  0  ~  90  degree  ) 

(Q)  Enter  the  azimuth  angle 

If  the  user  does  not  set  the  geophysical  parameters  and  specified  date  for  the 
tropospheric  model  by  passing  calling  arguments  to  I/O  subroutine  and  raytrace 
subroutine,  error  messages  will  be  displayed  and  will  prompt  input  questions  to  validate 
feature,  option,  and  other  previously  inputted  informations.  Figure  9  and  10  show  the 
syntax  for  the  interface  and  describe  calling  arguments. 

5.2.  Input  &  Output  Subroutine 

Figure  1 1  is  a  diagram  of  the  subroutine  which  makes  up  the  I/O  routine.  If  the 
user  selects  “ECM  or  HIRAS  or  MRF  database”  source  and  without  grid  cross  feature. 
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call  getwxfil  (GCrossFIg,  Source,  Lat,  Lon,  Press,  Temp,  ReIHum,  MsIFlg,  Mon,  Hour,  Day,  Yr,  LFCtropol, 
LFCtropo2,  Dataray,  Datarayflg) 


INPUTS: 

GCrossFIg 

(lnteger*2) 

0~Disable  Grid  Crossing  feature 
1-Enable  Grid  Crossing  feature 

Source 

(integer*2) 

Data  Source 

1- ECM 

2- MRF 

3- Hiras 

4- L)ser 

Lat 

(real*8) 

(Latitude  -90  to  +90) 

Lon 

(real*8) 

(Longitude  0  to  360  or -180  to  +180) 

^  Press 

(real*8) 

Surface  Pressure  (mbars) 

^  Temp 

(real*8) 

Surface  Temperature  (“IQ 

^  ReIHum 

(real*8) 

Surface  Relative  Humidity  (%) 

MsIFlg 

(integer*2) 

Reference 

0-Ellipsoid 

1-Mean  Sea  Level  (Default) 

Mon 

(integer*2) 

Month  (1  to  12) 

Hour 

(integer*2) 

Hour  1-OOOOHrs 

2- 0600Hrs 

3- 1200Hrs 

4- 1800Hr5 

Day 

(integer*2) 

Day  1  to  31 

Yr 

(integer*4) 

Year  1991, 1992,etc. 

LFCtropol 

(integer*2) 

Logical  File  Code  1  ( Unit  number ) 

LFCtropo2 

(integer*2) 

Logical  File  Code  2  ( Unit  Number ) 

OUTPUTS: 

Dataray 

(real*8) 

array(0;730,1:7) 

(grid) 

(lat) 

(Ion) 

(refractivity) 

(Height  Coefficient) 

(Range  Error  Variance) 

(Angie  of  Arrival  Error  Variance) 

DatarayFIg  (integer*2)  Error  flag 

O--N0  error  during  data  retrieval 
1 --Error  during  data  retrieval 


NOTE:  Superscript  ®  are  for  input  for  source  4  and  Output  for  data  source  ECMWF,  MRF,  and  HIRAS. 


Figure  9.  I/O  routine  syntax  and  calling  arguments 
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call  raytrace  (MsIFIg,  TmeZenAng,  TropoHgt,  StnAlt,  Az,  RangeSat,  c,  freq,  Dataray,  GCrossFIg,  STDFIg, 
TimeDelay,  TimeDelaySTD,  ZenLnchAng,  ZenLnchAngSTD,  RngDelay,  RngDelaySTD,  ErrFIg.DiffRange) 


INPUTS: 


MsIFIg 

(integer*2) 

Reference 

0--Ellipsoid 

1--Mean  Sea  Level  (Default) 

TrueZenAng 

(rears) 

Zenith  angle  to  line  of  site 

TropoHgt 

(rears) 

Height  of  Troposphere  (Default  27  km) 

StnAlt 

(rears) 

Station  Altitude  (m) 

Az 

(rears) 

Azimuth  angle(0  ~  360) 

RangeSat 

(real*8) 

Range  site  to  receiver  (Not  used  currently) 

c 

(rears) 

Speed  of  light  in  Vacuum 

freq 

(rears) 

Frequency  (Not  used  currently) 

Dataray 

(rears) 

array(0;730,1:7) 

(grid) 

(lat) 

(Ion) 

(refractivity) 

(Height  Coefficient) 

(Range  Error  Variance) 

(Angle  of  Arrival  Error  Variance) 


GCrossFIg 

(lnteger*2) 

0-Disable  Grid  Crossing  feature 

1 -Enable  Grid  Crossing  feature 

STDFIg 

(lnteger*2) 

0-No  Standard  Deviation  desired 

1-Standard  Deviation  desired 

OUTPUTS; 

TimeDelay 

(rears) 

Time  Delay  (ns) 

TimeDelaySTD 

(rears) 

Time  Delay  Variance  (ns^) 

ZenLnchAng 

(rears) 

Zenith  Launch  Angle  (Antenna  pointing  angle) 

ZenLnchAngSTD 

(rears) 

Zenith  Launch  Angle  Variance  (deg^) 

RngDelay 

(rears) 

Range  Error  (m) 

RngDelaySTD 

(rears) 

Range  Error  Variance  (m^) 

^rrFIg 

(integer*2) 

Error  flag 

0-No  error  during  raytrace 

1-Error  during  draytrace 

DiffRange 

(rears) 

Range  difference  exiting  troposphere(m) 

(Not  used  currently) 


NOTE;  Italic  Letter  is  not  used  currently 


Figure  10.  Raytrace  routine  syntax  and  calling  arguments 
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Figure  1 1  Diagram  of  I/O  subroutine 


II 
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Note:  Input  Parameter  Italic 

Output  Bold  Italic 
*  Input/Output 


this  routine  calls  grid  calculating  function.  If  the  returned  grid  from  the  grid  calculation 
function  is  out  of  range,  the  program  will  terminate.  If  not,  grid  number  will  go  to  the 
data  accessing  subroutine  to  extract  user’s  desired  data  with  respect  to  grid  number.  In 
order  to  test  whether  the  user  chooses  the  right  one,  I/O  subroutine  generates  print 
statement  for  verification  of  I/O  status.  When  opening  files  successfully  accesses  to  the 
desired  data,  the  lOSTAT  identity  in  the  read  statement  indicates  whether  the  file  has 
been  successfully  opened.  If  no  error  message  is  returned,  the  selected  data  is  stored  in  an 
array.  Thus,  I/O  routine  uses  dynamic  memory  allocation  method  based  on  specified 
geographic  location  for  data  storage.  Table  6  shows  dynamic  memory  allocation  size  for 
the  geophysical  location. 

A  ray  path  at  0°  elevation  angle  indicates  the  maximum  distance  from  ground 
station  to  the  target  or  the  spacecraft.  A  ray  path  distance  to  the  top  of  troposphere  (  27 
Km  from  the  surface)  is  approximately  587.49  Km  at  0°  elevation  angle.  In  table  6,  we 
can  see  the  diSerent  memory  size  based  on  latitude.  For  example,  the  user  selects  the 
latitude  range  between  60  ~  -60  degree,  I/O  program  allocates  array  size  (9,13).  If  the 
user  selects  “User  database”  source,  this  calls  user  info  subroutine  that  returns  surface 
pressure,  temperature,  and  humidity  needed  to  generate  surface  refractivity. 


Table  6.  Memory  allocation  for  different  latitude 


Latitude 

Memory  Size 

Note:  First  parameter  refers  to  latitude  for 
grid.  Second  refers  to  longitude  for  grid. 

80  ~  90  degree 

array  (  9,  73  ) 

70  ~  80  degree 

array  (  9,  25  ) 

60  ~  70  degree 

array  (  9,  17  ) 

60  ~  -60  degree 

array  (  9,  13  ) 

-60  ~  -70  degree 

array  ( 9, 17  ) 

-70  ~  -80  degree 

array  (  9, 25  ) 

-80  ~  -90  degree 

array  (  9,  73  ) 

If  inputs  are  invalid,  the  routine  will  be  run  out  of  range  error  message.  Given  that 
all  inputs  are  valid,  I/O  routine  reads  reference  coefficient  height  with  respect  to  input. 
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Each  subroutine  is  validated  by  using  example.  The  function  of  each  subroutine  is  as 
follows; 

a.  user  info 

This  subroutine  calls  cal  refrac  subroutine  that  calculates  surface  refractivity  and 
reads  reference  coefficient  height  data  file  containing  the  undulation  parameter  of  geoid 
and  reference  coefficients  height.  There  are  two  different  version  (  SUN  and  VAX  ).  In 
SUN  version,  the  operation  system  type  parameter(ostype)  is  equal  “sun”,  reads  the  data 
file  (  refht/refht.new)  under  user’s  installed  path.  In  VAX  version,  reads  the  data  file  ( 
[.refhtjrefht.vax  ).  In  order  to  validate  whether  reading  user’s  desired  data  file  or  not, 
added  printer  statement. 

Inputs  are  LFCtropol,  mslflg,  ostype,  date,  latitude,  longitude,  pressure, 
temperature,  and  relative  humidity.  Refi'activity  and  reference  height  are  stored  in  an 
array,  dataray.  Dataray  and  datarayflg  (error  flag)  are  passed  out.  Validation  input, 
expected  output  and  actual  output  from  subroutine  are  : 


Expected  Grid  Number 

Output  fi-om  subroutine 

23.5 

33.7 

1310 

1310 

-23.5 

-33.7 

6467 

6467 

33.6 

226.6 

1963 

1963 

b.  cal  refrac 

This  subroutine  calculates  surface  refractivity.  The  refractivity  for  air  which 
contains  water  vapor  is  given  as 


77.6  *  p  e  *  3.73  x  10^ 

N  =  -  •  +  - 

T  T^ 

e :  water  vapor 

p  :  atmosphere’s  barometric  pressure 
T :  atmosphere’s  temperature. 

The  first  part  of  the  right  hand  side  of  the  equation  computes  dry  refi'activity.  The  second 
part  computes  wet  refractivity.  Inputs  are  surface  pressure,  temperature,  and  relative 
humidity.  Output  is  the  total  surface  refractivity. 
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c.  direct  access 


This  subroutine  reads  data  files,  and  then  accesses  data  directly  to  datafiles.  Inputs 
are  cm,  source,  ostype,  grid,  date,  LFCtropol,  LFCtropo2.  Output  parameter  is  called  as 
parm  array.  In  order  to  validate  direct_access  subroutine,  inputs  are  passed  through 
user  info  subroutine  and  outputs  are  stored  in  parm  array.  When  the  input  and  output 
data  with  respect  to  user’s  input  do  not  match  ,  for  example,  the  input  grid  number  is  not 
equal  to  the  output  grid  number,  I/O  read  or  open  file  error  messages  are  displayed  and 
terminate  program. 

d.  grabGrid 

This  subroutine  organizes  an  array  of  data,  dataray,  in  an  organized  manner.  The 
organization  allows  for  faster  reading  of  data  files  by  creating  a  spread  sheet  that  stores 
grids  surrounding  initial  input  position  point.  Grids  are  stored  in  an  array  that  surrounds 
the  position  point  in  the  counterclockwise  direction.  Grids  spread  outward  from  initial 
point  in  rings  which  are  in  odd  square  matrix.  In  cases  where  initial  position  lies  between 
70  and  80  degrees  latitude,  the  grids  spread  in  odd  square  matrix  up  to  four  rings  then 
spreads  further  in  longitude  direction  adding  six  columns  of  grids  to  the  right  and  left  of 
the  square.  The  additional  columns  accounts  for  the  increased  possibility  of  grid  crossing. 
Inputs  are  cm,  date,  azimuth,  LFCtropol,  LFCtropo2,  mslflg,  latitude,  longitude,  source, 
ostype,  pressure,  temperature,  relative  humidity.  Outputs  are  dataray  and  datarayflg. 

Table  7  shows  the  result  of  validation  for  storing  data  in  an  array  to  different  user  input. 
The  array  size  are  relative  with  respect  to  table  6  in  section  5.2. 


Table  7  Dynamic  memory  allocation  of  different  user  input 


33  _  27 


m 

14 

13 

12 

11 

4 

3 

2 

10 

17 

5 

Tnpnt _ 

1 

9 

18 

6 

7 

8 

,  24 

19 

20 

21 

22 

23 

39 _ 45 


Note  ;  Each  number  represents  the  array  index. 
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5.3.  Raytrace  Subroutine 


There  are  three  main  parts  in  the  raytrace  subroutine.  Figure  12  shows  a  diagram 
of  the  raytrace  subroutine.  Part  one  is  an  iterative  approach  used  to  compute  the 
aggregate  refractive  effect  of  radio  propagation  from  ground  station  to  the  top  of 
troposphere  (  27  Km  from  the  surface).  Given  a  specific  refractivity-height  profile,  a  ray 
path  with  increasing  elevation  angles  are  ‘traced’  up  to  the  top  of  troposphere.  Geometric 
techniques  are  then  used  to  compute  the  true  elevation  angle  given  initial  LOS  angle.  An 
initial  estimate  of  the  actual  signal  elevation  angle  is  obtained  by  adding  an  angle  of  arrival 
error  which  was  computed  through  the  assumption  that  LOS  angle  is  an  initial  true 
elevation  angle.  Snell’s  law  for  spherical  geometry  together  with  the  refractivity  height 
profile  may  now  be  used  to  compute  a  distance  of  ray  path  at  each  height  in  the 
troposphere. 

The  basic  assumption  is  that  the  troposphere  is  considered  to  be  stratified  into  m 
height  profiles  and  constant  refractivity  profile  N™.  In  this  subroutine,  the  height  profiles 
is  created  by 

0  ~  1,000  meters  ;  increment  100  meters 

1,000  ~  10,000  meters  ;  increment  500  meters 

10,000  --  top  of  troposphere  :  increment  1,000  meters. 

Using  the  refractivity  and  the  height  profiles  reaching  up  to  the  tropospheric  height,  a  ray 
path  in  the  troposphere  can  be  obtained  using  Snell’s  law  and  Bouger’s  rule  for  spherical 
geometry  as  shown  in  the  following: 


sin  lo 


ro 


ao 


ri 


and 


no  ro  cos  ao  =  ni  ri  cos  ai 

where  ro  is  the  radius  of  the  Earth,  ri  =  ro  +  ho ,  n  is  the  radio  refractive  index  of  the 
atmosphere  at  any  layer,  and  a  is  the  apparent  elevation  angle  for  each  refractive  index 
profile,  and  I  is  an  incident  angle  between  different  two  layers.  The  ray  path  computed  at 
each  height  of  the  troposphere  can  be  used  to  compute  the  range  using  refractive  bending 
of  the  ray.  The  range  error  equation  is 


Figure  12  Diagram  of  Raytrace  subroutine 
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Note:  Input  Parameter  Italic 

Output  Bold  Italic 


where  the  distance  Rj  is  given  by 
Rf  =  r/  +  -  2  Tj  fj+i  cos  0j 

and  the  distance  Rom  is 


Rom  ro  Tm+l  ”  2  roTm+l  COS 


m 


H  0 

j=0 


and 

n 

0j  —  “ 

2 


The  refraction  angle  error  Aou  ,  which  is  the  difference  between  the  apparent  elevation 
angle  and  the  LOS  elevation  angle,  can  then  be  determined  from 


Aam  =  ao-aom. 


Part  one  of  the  program  uses  many  mathematical  equations  such  as  arccosine,  arcsine,  and 
square  root  function.  In  order  to  validate  above  mathematical  equations,  the  value  of 
arccosine  and  arcsine  is  assumed  to  be  greater  than  1.  Then,  it  returns  an  out-of-range 
error  message  and  resets  the  value  to  one.  If  the  square  root  error  occurs,  the  domain 
error  message  is  returned  and  program  will  be  terminated.  Part  two  of  the  program  uses 
climatology  weather  data  to  construct  refractivity  for  each  height  profiles,  defined  in  part 
one,  based  on  surface  data.  In  this  part,  it  calculates  each  comer  point  on  the  grid  which 
encloses  the  initial  location.  Then,  finds  the  shortest  perpendicular  distance  from  initial 
location  to  grid  lines.  It  calls  position  subroutine  to  find  the  last  projected  ground  location 
to  the  top  of  tropospheric  region.  If  the  grid  of  the  last  location  is  equal  to  the  initial  grid, 
raytrace  program  will  skip  directly  to  part  three  of  the  program.  If  not,  this  part  will 
check  grid  cross  index.  If  the  grid  cross  index  is  false,  this  part  will  be  skipped  and  will 
proceed  to  part  three.  If  the  grid  cross  index  is  true,  it  will  call  position  subroutine  to  find 
each  projected  location  and  grid  number  for  each  height  profiles.  If  the  grid  number  is  not 
equal  to  initial  grid,  it  will  access  surface  refractivity  for  different  grid  number  and 
compute  refractivity  with  respect  to  height  profile.  This  method  to  find  each  refractivity 
for  height  profile  will  be  iterated  until  height  reaches  the  troposheric  location.  Then 
refractivity  will  be  passed  into  the  third  part  of  the  program.  The  third  part  of  program 
traces  a  ray  from  the  Earth’s  ground  surface  to  the  top  of  troposphere  and  compute  range 
error,  time  delay,  and  angle  of  arrival  error. 
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a.  position 


This  subroutine  computes  coordinates  and  distances  for  curved  paths.  This  calls 
geodetic  subroutine  to  convert  given  initial  geodetic  location  to  earth-center  earth-fixed 
(ECEF)  coordination.  It  calculates  position  vector  (  south,  east,  and  zenith  component) 
fi'om  the  ground  station  to  the  satellite  point  using  true  elevation  and  azimuth.  Adding 
geodetic  vector  to  position  vector  returns  a  new  vector  position  over  the  Earth.  Then  the 
routine  calls  newwgs84latlon  subroutine  to  compute  the  subsatellite  point  on  Earth’s 
surface.  The  conversion  of  geodetic  to  geocentric  coordination  is  described  in  appendix 
A.  Inputs  are  gdlatl,  gdlonl,  altl,  azl,  ell,  hgol.  Outputs  are  new  geocentric  latitude 
and  longitude. 

b.  newwgs84latlon 

This  subroutine  computes  WGS-84  geodetic  latitude  and  longitude  from  ECEF 
coordinates.  It  computes  initial  geodetic  latitude  and  longitude  of  subsatellite  point. 
Then,  it  corrects  the  subsatellite  point  using  four  iteration  loop.  Input  is  ECF  position, 
xecf  Outputs  are  geodetic  subsatellite  latitude  “glaf’  degrees,  geodetic  subsatellite 
longitude  “glon”  degrees,  altitude  above  ellipsoid  “alt”,  and  distance  from  the  center  of 
Earth  to  the  surface  point  “rmag”. 

c.  geodetic 

This  subroutine  computes  Earth-Centered  Earth-Fixed  (ECEF)  cartesian 
coordinates  given  geodetic  latitude,  longitude  and  altitude.  Inputs  are  geodetic  latitude 
“xlat”  in  degrees,  longitude  “xlon”  in  degrees,  and  height  “xh”  above  the  reference 
ellipsoid  in  kilometers.  Outputs  are  ECF  x-coordinate  “xx”,  y-coordinate  “yy”,  z- 
coordinate  “zz”  in  kilometers. 

6.  CONCLUSION 

Four  databases  are  selected  as  most  valuable  sources  for  climatological  data.  The 
7  -  variable  surface  data,  is  derived  fi'om  the  17  variable  multilayer  data  for  use  for  the 
Modified  Exponential  Model.  This  model  builds  a  multilayer  atmosphere.  The 
constructed  atmospheric  data  are  then  used  with  Millman’s  Stratified  layer  method  to 
calculate  range  error,  time  delay,  and  angle  of  arrival  error.  Tests  performed  on  130 
various  worldwide  location  indicate  that  most  realistic  results  are  obtained  using  the 
Modified  Exponential  Model.  Numerous  models  are  examined  including  Hopfield(  with 
surface  weather  conditions  )  and  ESS  models.  Error  parameter  statistics  show  that  errors 
are  less  than  six  percent  standard  deviation  in  time  delay  and  range  error,  with  less  than  24 
percent  standard  deviation  in  zenith  launch  angle.  In  addition,  azimuth  sensitivity  tends  to 
be  a  function  of  location  and  time. 

The  Modified  Exponential  Model  source  code  is  created  to  be  a  user-fiiendly  for 
the  needs  of  radio  propagation  system  engineer  who  need  many  propagation  phenomena 
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analysis.  Data  files  are  provided  in  ASCII-format  file  and  can  be  easily  reformatted  for 
direct-access.  Although  great  care  should  be  taken  in  modifying  the  program,  the  author 
would  appreciate  any  comment  regarding  for  improving  the  software  quality. 
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Appendix  A 

Conversion  Geodetic  to  Earth-Center  Earth-Fixed  (  ECEF) 
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There  are  two  most  commonly  used  definitions  of  latitude.  Figure  13  shows  Geocentric 
and  Geodetic  latitude. 


ae  =  Equatorial  radius  ( 6378. 167  Km ) 
be  =  Polar  radius  (6356.785  Km) 


Figure  13  Geocentric  and  geodetic  latitude 


The  angle  L’  is  called  “geocentric  latitude”  and  is  defined  as  the  angle  between  the 
equatorial  plan  and  the  radius  fi-om  the  geocenter.  The  angle  L  is  called  “  geodetic 
latitude”  and  is  defined  as  angle  between  the  equatorial  plan  and  the  normal  to  the  surface 
of  the  ellipsoid.  The  word  “latitude”  usually  means  geodetic  latitude. 

What  we  need  now  is  a  method  of  calculating  the  station  coordinates  of  a  point  on  the 
surface  of  our  reference  ellipsoid  when  we  know  the  geodetic  latitude  and  longitude  of  the 
point  and  its  height  above  mean  sea  level.  Consider  an  ellipse  comprising  and  rectangular 
coordinate  system  as  shown  in  Figure  14  . 
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Figure  14  Station  Coordinates 


We  will  first  determine  the  X  and  Z  coordinates  of  a  point  on  the  ellipse  assuming  that  we 
know  the  geodetic  latitude  L.  It  will  then  be  a  simple  matter  to  adjust  these  coordinates 
for  a  point  which  is  a  known  elevation  angle  above  the  surface  of  the  ellipsoid  in  the 
direction  of  the  normal.  It  is  convenient  to  introduce  the  angle  P  ,  the  “reduced  latitude” 
which  is  illustrated  in  Figure  14  .  The  X  and  Z  coordinates  can  immediately  be  written  in 
term  of  P  if  we  note  that  the  ratio  of  the  z  coordinate  of  the  point  on  the  ellipse  to  the 
corresponding  Z  ordinate  of  a  point  on  the  circumscribed  circle  is  just  hja^.  Thus, 


be 

z  = - ae  sin  P 

ae 

But,  for  any  ellipse ,  a^  =  and  e  =  c/a,  where  e  is  eccentricity,  so 
be  =  ae  Vl-e^ 
and 

z  =  ae  Vl-e^  sin  3 


We  must  now  express  sin  P  in  terms  of  the  geodetic  latitude  L  and  the  constants  ae  and  be. 
From  elementary  calculus  we  know  that  the  slope  of  the  tangent  to  the  ellipse  is  just  dz/dx 
and  the  slope  of  the  normal  is  -dx/dz.  Since  the  slope  of  the  normal  is  just  tan  L,  we  can 
write 


dx 

tan  L=  -  - 

dz 

The  differentials  dx  and  dz  can  be  obtained  by  differentiating  the  expressions  for  x  and  z 
above.  Thus, 

dx  =  -  a*  sin  P  dp 

dz  =  a«  Vl-e^  cos  P  dp 

and 


tanL  = 


tan  3 


or 


tanP  =  Vl-e^  tanL  = 


Vl-e^  sin  L 


cosL 
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Suppose  we  consider  this  last  expression  as  the  quotient 


ae  (  1  -  e^  )  sin  L 
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Appendix  B 

Interface  syntax  for  each  subroutine 
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Subroutine:  DIRECT  ACCESS 


Calling  Syntax: 

CALL  DIRECT_ACCESS(TDATE,  SOURCE,  OSTYPE,  GRID,  DATE, 
LFCtropol,LFCtropo2,  FARM,  REFHT,  UNDUE) 

Description  of  Arguments 

Arguments  passed  from  getwxfrl  subroutine  called  in  Main: 


TDATE 

SOURCE 

OSTYPE 

GRID 

DATE 

LFCtropol 

LFCtropo2 


[character*70] 

Data  source(ecm,mrf,hiras,user)  [integer*2] 

System  type  (sun/vax)  [character*  3] 

Grid  number  [integer*4] 

Date(4)  month,hour,day,year  in  an  array  [integer*2] 
Logical  File  Code  1  [integer*2] 

Logical  File  Code  2  [integer*2] 


Arguments  returned  to  getwxfrl  subroutine: 


FARM  :  Parm(  1:6),  parameter  containing  information  relative 

to  the  grid  i.e.  pressure,  temperature,  etc  [real*8] 
REFHT  :  Reference  height  [real*  8] 

UNDUE  :  undulation  [real*8] 
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Subroutine;  GEODETIC 


Calling  Syntax: 

CALL  GEODETIC(XLAT,  XLON,  XH,  XX,  YY,  ZZ,  RR) 

Descriiption  of  Arguments 

Arguments  passed  from  position  subroutine  called  in  Main; 

XLAT  :  Geodetic  latitude  in  degree  [real*  8] 

XLON  :  Geodetic  longitude  in  degree  [real*8] 

XH  :  Height  above  reference  ellpisoid  in  kilometer  [real*8] 

Arguments  returned  to  position  subroutine: 

X  ECF  coordinate  in  kilometers  [real*  8] 

Y  ECF  coordinate  in  kilometers  [real  *8] 

Z  ECF  coordinate  in  kilometers  [real*  8] 

Radial  distance  from  earth  center  to 
(x,y,z)  position  in  kilometer  [real*8] 


XX 

YY 

ZZ 

RR 
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Subroutine:  GRABGRID 


Calling  Syntax: 

CALL  GRABGRID(TDATE,  DATE,  AZ,  LFCtropol,  LFCtropo2, 
MslFlg,  LAT,  LON,  SOURCE,OSTYPE,  PRESS,  TEMP,  RELHUM, 
DATARAY,  DATARAYFLG) 

Descriiption  of  Arguments 

Arguments  passed  from  getwxfil  subroutine  called  in  Main: 


TDATE 

DATE 

AZ 

LFCtropol 

LFCtropo2 

MslFlg 

LAT 

LON 

SOURCE 

OSTYPE 

PRESSURE 

TEMP 

RELHUM 


:  [character*70] 

:  Date(4)  month,hour,day,year  in  an  array  [integer*2] 
:  Azimuth  angle(0-360)  [real*  8] 

:  Logical  File  Code  1  [integer*2] 

:  Logical  File  Code  2  [integer*2] 

:  Reference  (ellipsoid/mean  sea  level)  [integer*2] 

:  Latitude  in  degree  [real  *8] 

:  Longitude  in  degree  [real*8] 

:  Data  source(ecm,mrf,hiras,user)  [integer*2] 

:  System  type  (sun/vax)  [character*3] 

:  Surface  pressure  in  millibar  [real*  8] 

:  Surface  temperature  in  Kelvin  [real*  8] 

:  Surface  relative  humidity  [real*8] 


Arguments  returned  to  getwxfil  subroutine: 


DATARAY  :  Dataray(0:730, 1:8),  an  array  [real*8] 

Contains  following: 

grid  number,  latitude,  longitude,  refractivity, 
height  coeflBcient,  range  error  variance, 
angle  error  variance 
DATARAYFLG  :  Error  flag 

0--No  error  during  data  retrieval 
1 --error  during  data  retieval 
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Subroutine:  NEWWGS84LATLON 


Calling  Syntax; 

CALL  NEWWGS84LATLON(XECF,  NEWLAT,  NEWLON,  ALT,  RMAG) 
Descriiption  of  Arguments 


Arguments  passed  from  position  subroutine  called  in  Main; 
XECF  :  XECT(7),  ECF  position  [real*8] 


Arguments  returned  to  position  subroutine; 


NEWLAT 

NEWLON 

ALT 

RMAG 


Geodetic  subsatellite  latitude  in  degree  [real*8] 
Geodetic  subsatellite  longitude  in  degree  [real*8] 
Altitude  above  ellopsoid  [real*8] 

Distance  from  center  to  subsatellite  point  [real*8] 
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Subroutine:  POSITION 


Calling  Syntax: 

CALL  POSITION(gdLATl,  gdLONl,  ALTl,  AZl,  ELI,  RHOl, 
NEWLAT,  NEWLON) 

Descriiption  of  Arguments 


Arguments  passed  from  raytrace  subroutine  called  in  Main: 


gdLATl 

gdLONl 

ALTl 

AZl 

ELI 

RHOl 


:  Input  latitude  in  degree  [real*  8] 

:  Input  longitude  in  degree  [real*8] 

:  Input  altitude  in  meter  [real*8] 

:  Input  azimuth  in  degree  [real*8] 

:  Input  elevation  in  degree  [real*8] 

:  Input  ray  path  in  kilometers  [real*8] 


Arguments  returned  to  raytrace  subroutine: 


NEWLAT  :  New  latitude  position  in  degree  [real*  8] 

NEWLON  :  New  longitude  position  in  degree  [real*8] 

Errflg  :  Error  flag 
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Subroutine;  USER  INFO 


Calling  Syntax: 

CALL  USER_INFO(DATE,  LFCtropol,  MslFlg,  LAT,LON,  PRESS, 
TEMP,  RELHUM,  OSTYPE,  DATARAY,  DATARAYFLG) 

Description  of  Arguments 

Arguments  passed  from  getwxfil  subroutine  called  in  Main; 


DATE 

LFCtropol 

MslFIg 

LAT 

LON 

PRESS 

TEMP 

RELHUM 

OSTYPE 


Date(4)  month,  hour,day,  and  year  in  array  [integer*!] 
Logical  File  Code  1  [integer*!] 

Reference  (ellipsoid/mean  sea  level)  [integer*!] 
Latitude  in  degree  [real*8] 

Longitude  in  degree  [real*8] 

Surface  pressure  in  millibar  [real*8] 

Surface  temperature  in  Kelvin  [real*8] 

Surface  relative  humidity  in  percent  [real  *8] 

System  type  (sun/vax)  [character*3] 


Arguments  returned  to  getwxfil  subroutine: 


DATARAY  :  Dataray(0; 730,1; 8),  an  array  [real*8] 

Contains  following; 

grid  number,  latitude,  longitude,  refractivity, 
height  coefficient,  range  error  variance, 
angle  error  variance 
DATARAYFLG  ;  Error  flag 

0~No  error  during  data  retrieval 
1— error  during  data  retieval 
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Appendix  C 


Flow  Charts  of  Data  Processing  and  Source  Code 
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Data  Processing  Flow  Diagram  for  Tropospheric  Propagation  Database 
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Real  Time 
User  Input  Data 


Flow  Chart  of  Main  Program 


jS> 
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End 


Flow  Chart  of  Input/Outout  Program 
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Flow  Chart  of  Raytrace  Program 
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Time  delay.  Range  error,  and  Angle  of  arrival  error 


Appendix  D 

Percent  Standard  Deviations  Time  Delay 
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ECMWF  Time  Delay  Percent  Standard  Deviation  Contour 


o 
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ECMWF  Time  Delay  Percent  Standard  Deviation  Contour  Map  for  Jul  (Angle  =  0  deg) 
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Appendix  E 


Azimuth  Sensitivity  HlRAS  data 
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HIRAS  Percent  Change  in  Time  Delay  (Feb  0600  Surface  Data,  45  deg  Azimuth  Increments,  Angle  =  0  deg) 
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HIRAS  Percent  Change  in  Time  Delay  (Aug  0600  Surface  Data,  45  deg  Azimuth  Increments,  Angle  =  0  deg) 
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Percent  Change  in  Time  Delay  (Nov  0600  Surface  Data,  45  deq  Azimuth  Increments 
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